#include "func.h"
#include "is_eq.h"
#include "int2str.h"
#include "zadacha.h"
#include "length.h"
#include "rk4.h"
#include "rkf45.h"
#include "table.h"

#include <fstream>
#include <iostream>
#include <iomanip>
using namespace std;


double func (std::valarray<double> w, file_names& Fnames, input& data) {
    unsigned m = data.q,    //размерность системы: x + z
             iter_i = 0,
             mnoj = 100,
             count = 0,
             N = data.N;
    double coeff0 = 0.9, coeffu = -0.2, coeffl = -0.25,
           delta = 0,
           max_oshibka = 0,
           t = data.t,
           tf = data.tf,
           th = fabs(tf - t)/N;

    unsigned max_iter = N*mnoj;

    std::valarray<double> x(m);
    x = data.x;                 //значения на левом конце


    valarray<double> x1(m);

    std::vector<double> tabl;   //итоговые результаты
    std::ofstream out(Fnames.runge);
    tabl.push_back(t);
    for (unsigned i=0; i<m; i++) tabl.push_back(x[i]);

    switch (data.tip_integr) {
        case 1:
              for (unsigned i=0; i<N; i++) {
                  rk (t, th, &x[0], w, data);
                  t += th;
                  tabl.push_back(t);
                  for (unsigned i=0; i<m; i++) tabl.push_back(x[i]);
              }
              break;
        case 2:
              while (t < tf) {
                  count++;
                  if (count>max_iter) {
                      cout << "Превышено максимальное число итераций" << endl;
                      return I(x,data);
                  }

                  max_oshibka = rkf45 (t, th, x, w, &x1[0], data);
                  //cout << "max_oshibka= " << max_oshibka << endl;
                  if (max_oshibka<=1) {
                      t += th;
                      delta = coeff0*pow(max_oshibka, coeffu);
                      if (delta>4) {
                          th *= 4;
                      }
                      else if (delta>1) th *= delta;
                      if (t+th > tf) th = tf - t;
                      x = x1;
                      //cout << "th= " << th << endl;
                      tabl.push_back(t);
                      for (unsigned i=0; i<m; i++) tabl.push_back(x[i]);
                  }
                  else {
                      delta = coeff0*pow(max_oshibka, coeffl);
                      if (delta < 0.1) th *= 0.1;
                      else th *= delta;
                  }
              }
              iter_i = count;
              break;
    }

    //вывод траекторий: t x
    for (unsigned j=0; j<tabl.size(); j++) {
        out << std::setprecision(10)
            << std::setiosflags(std::ios::fixed | std::ios::showpoint)
            << std::setw(15) << std::setiosflags(std::ios::left) << tabl[j];
        if ((j+1)%(m+1)==0) out << std::endl;
    }

    for (unsigned i=0; i<data.m; i++) {
        out << " otklon[" << i << "]= " << data.df[i];
    }
    out << endl;

    out.close();

    return I(x,data);
}

std::valarray<double> upr (double t, std::valarray<double> w, input& data) {
    unsigned num_upr = data.num_upr;
    std::valarray<double> c(num_upr);

    c = data.c;

    std::valarray<std::valarray<double> > tau;
    tau.resize(num_upr);

    for (unsigned i=0; i<num_upr; i++) tau[i].resize(c[i]+2);
    tau = table(w, data);

    std::vector<unsigned> sum(num_upr+1);
    sum[0]=0;
    const unsigned p=2;
    std::valarray<double> res_u(num_upr);


    unsigned n = data.n;//w.size();
    switch (data.tip_metoda) {
        case 1:
              for (unsigned i=1; i<num_upr+1; i++) {
                  sum[i] = (c[i-1]+1)+sum[i-1];
              }
              for (unsigned u_ind=0; u_ind<num_upr; u_ind++) {
                  for (unsigned i=0; i<c[u_ind]+1; i++) {
                      if (((t>tau[u_ind][i])||IsEq(t,tau[u_ind][i],1e-11)) && ((t<tau[u_ind][i+1])||(IsEq(t,tau[u_ind][i+1],1e-11)))) {
                          res_u[u_ind] = w[i+sum[u_ind]];
                      }
                  }
              }
              break;
        case 2:
              for (unsigned i=1; i<num_upr+1; i++) {
                  sum[i] = (c[i-1]+2)+sum[i-1];
              }
              for (unsigned u_ind=0; u_ind<num_upr; u_ind++) {
                  for (unsigned i=0; i<c[u_ind]+1; i++) {
                      if (((t>tau[u_ind][i])||IsEq(t,tau[u_ind][i],1e-11)) && ((t<tau[u_ind][i+1])||(IsEq(t,tau[u_ind][i+1],1e-11)))) {
                          res_u[u_ind] = (((w[i+1+sum[u_ind]]-w[i+sum[u_ind]])*t-(tau[u_ind][i]*w[i+1+sum[u_ind]]-tau[u_ind][i+1]*w[i+sum[u_ind]]))/(tau[u_ind][i+1]-tau[u_ind][i]));
                      }
                  }
              }
              break;
        case 3:
              for (unsigned i=1; i<num_upr+1; i++) {
                  sum[i] = (c[i-1]+1)+sum[i-1];
              }
              for (unsigned u_ind=0; u_ind<num_upr; u_ind++) {
                  for (unsigned j=0,i=0; (j<c[u_ind]+1) && (i<n); j++,i+=p) {
                      if (((t>tau[u_ind][j])||IsEq(t,tau[u_ind][j],1e-11)) && ((t<tau[u_ind][j+1])||(IsEq(t,tau[u_ind][j+1],1e-11)))) {
                          res_u[u_ind] = (((w[i+1+2*sum[u_ind]]-w[i+2*sum[u_ind]])*t-(tau[u_ind][j]*w[i+1+2*sum[u_ind]]-tau[u_ind][j+1]*w[i+2*sum[u_ind]]))/(tau[u_ind][j+1]-tau[u_ind][j]));
                      }
                  }
              }
              break;
    }
    return res_u;
}
